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VLTI/MIDI 10 MICRON INTERFEROMETRY OF THE FORMING MASSIVE STAR W33A 
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ABSTRACT 

We report on resolved interferometric observations with VLTI/MIDI of the massive young stellar 
object (MYSO) W33A. The MIDI observations deliver spectrally dispersed visibilities with values 
between 0.03 and 0.06, for a baseline of 45m over the wavelength range 8-13/j,m. The visibilities 
indicate that W33A has a FWHM size of approximately 120 AU (0.030") at 8/im which increases to 
240 AU at 13/im, scales previously unexplored among MYSOs. This observed trend is consistent with 
the temperature falling off with distance. ID dust radiative transfer models are simultaneously fit to 
the visibility spectrum, the strong silicate feature and the shape of the mid infrared spectral energy 
distribution (SED). For any powerlaw density distribution, we find that the sizes (as implied by the 
visibilities) and the stellar luminosity are incompatible. A reduction to a third of W33A's previously 
adopted luminosity is required to match the visibilities; such a reduction is consistent with new high 
resolution 70/im data from Spitzer's MIPSGAL survey. We obtain best fits for models with shallow 



dust density distributions of r and 



,-1.0 



and for increased optical depth in the silicate feature 



produced by decreasing the ISM ratio of graphite to silicates and using optical grain properties by 
Ossenkopf et al. (1992). 

Subject headings: stars: formation — stars: early type — techniques: interferometric 



1. INTRODUCTION 

Massive young stellar objects are luminous, embed- 
ded infrared (IR) sources that show many signs that 
they are still actively accreting mass. Their luminos- 
ity {L > lO'^L©) is such that they are expected to be 
ionizing their surroundings to produce an H II region, 
yet they only have weak radio emission due to ionized 
winds or jets (Hoare 2002). Most appear to be driving 
bipolar molecular flows and (sub-)millimetre interferom- 
etry is beginning to reveal evidence of rotating, disc-like 
structures on scales of hundreds of AU (Patel et al. 2005; 
Beltran et al. 2006; Torrelles et al. 2007). There is great 
interest in knowing the distribution of the infalling and 
outflowing material on smaller scales to provide clues 
to which physical processes are controlling the dynam- 
ics and setting the final mass of the star (Beuther et al. 
2007; Hoare et al. 2007). 

The mid-IR (8-13/im) emission from MYSOs is thought 
to arise in the warm (^--^300 K) dust in the envelope heated 
by the young star. Previous modelling (e.g. Churchwell 
et al. 1990) has indicated that the size of the mid-IR 
emission region should be unresolved at the typical dis- 
tances of MYSOs by single-dishes and this is borne out 
by observations (Kraemer et al. 2001; Mottram et al. 
2007). Exceptions occur when MYSOs are viewed close 
to edge-on when a dense torus can completely obscure 
the bright central region and only the warm dust in the 
outflow cavities is seen (de Buizer 2006). 

Here we present new results from mid-IR interferomet- 
ric observations of the MYSO W33A using the VLTI 
MIDI instrument. W33A has a kinematic distance of 
3.8 kpc and a luminosity as derived from IRAS fluxes 
of 1 X 1O^L0 (Faundez et al. 2004). It has weak, com- 
pact, optically thick radio continuum emission (Rengara- 
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jan & Ho 1996; van der Tak & Menten 2005) and broad 
('^lOO km s~^), single-peaked H I recombination emis- 
sion lines (Bunn et al. 1995) consistent with an ionized 
stellar wind origin. IR images from 2MASS and Spitzer's 
GLIMPSE survey clearly show a large scale monopolar 
nebula emerging to the SE, which is most likely to de- 
marcate the blue-shifted lobe of a bipolar outflow. 

The cold dust envelope has been studied at a few arc- 
second resolution at millimetre wavelengths (van der Tak 
et al. 2000), probing the imprint of the star formation 
process on the circumstellar dust distribution. Many dif- 
ferent scenarios have been put forward predicting this 
distribution. In case of W33A, Giirtler et al. (1991) mod- 
elled the SED with a constant density spherical envelope 
and found an inner cavity radius of 135 AU (35 mas at 
3.8 kpc) with the dust close to the sublimation tempera- 
ture. MIDI provides data on comparable size scales of 50 
mas, and since the mid-IR size will be somewhat larger 
than the sublimation zone we expect the emission to be 
well-resolved. 

2. OBSERVATIONS AND DATA REDUCTION 

MIDI is the VLTFs two telescope beam combiner that 
operates in the thermal IR (Leinert et al. 2003). W33A 
was observed with MIDI on 16 September 2005, employ- 
ing the UT2-UT3 telescope configuration for a projected 
baseline of 45.5 m and position angle of 47° east of north; 
this is perpendicular to the star's outflow direction. An 
interferometric calibrator star HD 169916 was observed 
during the same night. Lunar occultation and indirect 
methods reveal the calibrator diameter to be a few mas 
(Nather 1972; Pasinetti-Fracassini 2001), too small to be 
resolved by MIDI. Here we present visibilities between 
8-13/xm that are spectrally dispersed at a resolution of 
R = 30. The observations were executed in the High- 
Sens MIDI mode (for details see Przygodda et al. 2003; 
Chesneau et al. 2005; Leinert et al. 2004). We summa- 
rize the main elements of the procedure that produces the 



Fig. 1. — Left: MIDI (lower full line) and ISO flux calibrated spectra of W33A. Right: The correlated flux determined with EWS (dots 
with errorbars), and MIA (full line). Note the ragged behaviour of MIA at low flux levels. EWS errorbars correspond to the variations in 
visibilities as function of reduction parameters. EWS and MIA are consistent except at correlated flux levels below O.lJy (dashed line), 
these values must be considered upper limits (see Jaffe et al. 2004; Matsuura et al. 2006). 



raw dataset. The two beams are interfered producing 
two complementary interferometric channels that have 
a phase difference of n radians. The dispersed fringes 
are modulated according to an introduced optical path 
delay (OPD). Intcrfcrograms are recorded for a range 
in OPDs, a few millimetre around optical path length 
equalization, called a fringe scan. For both W33A and 
HD 169916 a total of 8000 interferograms, corresponding 
to 200 scans were recorded. Subtracting the two interfer- 
ometric beams eliminates the background and enhances 
the fringe signal. The coherent flux was extracted using 
two different MIDI software reduction packages: EWS 
(Jaffe 2004) and MIA (Kohler et al. 2005). MIA esti- 
mates the amplitudes in the power spectra of the fringe 
signal Fourier transform (incoherent estimation). The 
amplitudes are proportional to the correlated flux. EWS 
on the other hand aims at adding the fringes to maximize 
the signal to noise (coherent estimation). This method 
corrects first the fringe spectra for their corrupted phase 
caused by atmospheric and instrumental effects. EWS 
uses the fringe scan as phase reference to estimate the 
group delay due to the atmosphere. The observations 
show that the atmospheric group delay varied within a 
range of 50/im over 100s, indicating relatively good at- 
mospheric circumstances during the night. The Fourier 
transform of the group delay function reveals the typical 
sawtooth phase change due to the introduced OPD, in- 
dicating that fringes have been measured. Removing the 
atmospheric and the (known) instrumental group delays 
constitutes a linear correction to the dispersed fringe sig- 
nal, and straightens the dispersed fringe spectra, i.e. the 
phase is independent of wavelength. Next, the phase off- 
set due to varying water refraction index between the 
time of recording of the fringe spectra has to be ac- 
counted for. In principle all spectra can then be added to 
a final fringe spectrum. Final visibilities are obtained by 
taking two spectra of the source (one with each telescope) 
immediately after the interferometric measurement. The 
accuracy of the final visibility spectrum is limited by the 
sky brightness variation between the interferometric and 
flux measurement, and can amount to 10-15%. 

The flux spectra are corrected for sky contribution us- 
ing a median sky subtraction. Spectra are extracted by 
summing the counts within Str of the mean position, de- 
termined from a Gaussian fit to each column along the 



wavelength axis. This procedure results in a W33A spec- 
trum with an S/N - 100 and S/N ~ 300 for the stan- 
dard star. Four spectra are recorded, 2 for each telescope 
beam and combined via a geometric mean. This quantity 
is also what is obtained for the correlated fiux after beam 
combination and thus ensures consistency when deriving 
visibilities. The absolute flux calibration is done using 
HD 169916, which has an average flux density of 30 Jy 
between 8/im and 13/im (Cohen et al. 1999). The dif- 
ference in airmass between the observation of the flux 
calibrator and W33A is 0.05 leading to a negligible cor- 
rection on the flnal fluxes. Errorbars on the spectrum 
in Fig. [1] indicate the systematic difference in flux levels 
between the two telescopes beams. Absolute flux cali- 
bration is uncertain upto at least 35% due to difference 
in flux level of HD 169916 in each telescope beam. 

3. RESULTS 

Fig. [T] presents the MIDI flux spectrum and correlated 
flux spectrum from EWS and MIA. The flux spectrum 
is compared to the ISO-SWS spectrum, taken from the 
ISO Data Archive maintained by ESA. Both spectra are 
dominated by a very strong silicate feature, that con- 
tains solid-state ammonia and methanol absorption fea- 
tures, at 9 and 9.7/im respectively (see Gibb et al. 2000). 
At the central wavelength of the feature no flux was 
recorded, and the actual depth is unknown. The MIDI 
and ISO spectra portray a similar overall shape but the 
MIDI spectrum has a flux level that is about a factor 2 
less. In addition to the uncertainties due to flux calibra- 
tion, we ascribe this difference to the much larger ISO 
beam (80") in comparison to MIDI (2"). 

The corresponding visibility spectrum is obtained by 
dividing the correlated flux by the flux spectrum and is 
presented in left panel of Fig. [51 Visibilities are not plot- 
ted when corresponding to correlated fluxes smaller than 
O.lJy, a value below which the measurements become un- 
rehable (Jaffe et al. 2004, Matsuura et al. 2006). The 
visibility spectrum shows that the silicate wings are not 
strongly affected by the decrease in flux but instead fol- 
low the declining trend of the continuum visibilities. If 
we would represent the emission by a Gaussian emitting 
distribution then the FWHM size increases from 30 mas 
at 8/im to 60 mas at 13/im. 

4. RADIATIVE TRANSFER MODELLING 
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Fig. 2. — DUSTY model fits to visibilities and SED. Model is described by a r " density law, and Ay = 100, star parameters are 
L = 410'*Lq, T^ff = lOOOOK, and a distance of 3800 pc. Left: Calibrated MIDI visibilities of W33A (dots with errorbars). The level of 
systematic uncertainty is indicated by the bar in the right top corner. Black squares connected by the dashed line correspond to the best-fit 
model. Right: SED of W33A traced by ISO-SWS and LWS spectra and IRAS points (open circles). Full line is the model fit. IRAS/ISO 
large beam effect at the longer wavelength is illustrated by the improved photometry using MIPSGAL images (filled circle). 



Previous model fits of MYSO SEDs have indicated that 
simple spherical radiative transfer models with roughly 
constant densities best fit the near-IR and far-IR data 
(Churchwell et al. 1990; Giirtler et al. 1991). How- 
ever, evidence exists that the outer envelopes (10,000 AU 
scales) have steeper density profiles (Hoare et al. 1990; 
van der Tak et al. 2000). We explore here whether 
the visibilities produced by material on 100 AU scales (a 
scale previously unexplored), and the SED can be simul- 
taneously matched by simple ID spherically symmetric 
dust models. Arguably, such models are inadequate for 
MYSOs which are likely to consist of circumstellar disks, 
bipolar cavities, etc. Despite this, we explore whether 
the basic levels and trends of the dispersed visibilities 
can be matched by a single unresolved star deeply em- 
bedded in a dusty envelope before introducing more free 
parameters. 

For this purpose, we employ DUSTY, a code that 
solves the scaled ID dust radiative transfer problem (see 
Ivezic & Elitzur 1997). We used a spherically symmet- 
ric dust distribution illuminated by a central, unresolved 
star. The only non-scaled parameter entering the code 
is the dust sublimation temperature. Gas emission is 
not taken into account by DUSTY. Solutions are inde- 
pendent of the central source luminosity, and the SED 
and visibilities can be scaled accordingly. The luminos- 
ity is the prime stellar parameter that sets the inner dust 
sublimation radius and thus the size scale; an increase 
causes the size of the emitting region to be larger, as 
Tsub oc ^/L. The canonical parameters for W33A are 
L = 1 IO^Lq, 3.8 kpc (Faundez et al. 2004) and we ini- 
tially chose Toff = 25000 K typical for 15 Mq star. These 
values make the underlying star larger than a main se- 
quence Bl star, consistent with the idea that a star un- 
dergoing accretion at a high rate is swollen (e.g. Behrend 
& Maeder 2001). We tuned the model output to the ob- 
servations by varying DUSTY's input parameters, viz. 
radial density distribution, dust properties, and source 
Teff. The challenge is to construct a single model that 
fits the observed visibilities, silicate wings, and mid-IR 
SED up to lOO^m. We restrict the SED fitting to this 
wavelength interval as it is well-known that ID models 
have particular difhculty in reproducing the short wave- 
length region because of the sensitivity to the viewing 



angle (cf. Yorke & Sonnhalter 2002). 

We adopt a dust sublimation temperature of 1500 K 
and an MRN-DL dust mixture (Mathis et al. 1977; 
Draine & Lee 1984) with a typical interstellar graphite- 
silicate composition. The outer bound of the model is 
set at 1000 times the dust sublimation radius, where the 
temperature corresponds to the presumed ambient tem- 
perature of the ISM, between 10 and 25 K. The precise 
value for the outer radius does not affect our conclusions 
here. 

Previous studies have shown that W33A is best de- 
scribed by an Ay between 100 and 200 (Capps et al. 
1978; Giirtler et al. 1991). A first result is that the MIDI 
visibilities cannot be reproduced for any density distri- 
butions with power exponents between —2.0 and 0.0 for 
Av between 100 and 200 and L = 1 10^ Lq. The model 
visibilities at 10/zm are too small by an order of magni- 
tude. Decreasing Ay is not a solution for normal dust, 
because of W33A's exceptionally strong silicate absorp- 
tion feature. Average model sizes at lO/im of ~ 200AU 
are reached only if the luminosity is reduced to about 
a third. Such a reduction is justified when we consider 
70/zm MIPSGAL data that we have obtained from the 
Spitzer archive. These data are taken at a vastly superior 
resolution compared to both IRAS and ISO and reveal 
the presence of at least three point sources and strong 
diffuse emission within the IRAS beam. For W33A, we 
measure a 70/im flux of 1.1 10'^ Jy (20% uncertainty), as 
shown in Fig.O This indicates that the true luminosity 
is significantly less than deduced from IRAS data. 

Even with the luminosity reduced to 410'*Lo, ther- 
mally supported cores with dependency (Larson 
1969) are incompatible with the observed trend of de- 
creasing visibilities with wavelength and the SED. In- 
compatible models are also found if we adopt a constant 
infall velocity type distribution with r^^-^ (Shu 1977). 
Reasonable fits to the SED that produce sizes compa- 
rable to the MIDI visibilities are found for r"^-^ loga- 
tropic distributions (e.g. McLaughlin & Pudritz 1996), 
but again they do not reproduce the observed decreas- 
ing visibility trend with wavelength. DUSTY produces 
smaller sizes, and thus larger visibilities, for longer wave- 
lengths. This affects especially the red wing of the sili- 
cate feature, due to the diminishing opacity further out 
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in the silicate wing. This mismatch in visibihties is how- 
ever removed if we reduce the slope to a much shallower 
density distribution with a r~°'^ law or a constant den- 
sity law. These distributions fit the short wavelength 
region of the SED increasingly worse because of the loss 
of warm dust. These shallower density distributions also 
require the luminosity of the central object to be fainter 
than sB9observed. We thus find that for standard ISM 
dust, the SED and the sizes at scales of 100 AU limit the 
dust to follow distributions between r"^ *^ or r"*^'^ pow- 
erlaws. The required optical depths are relatively high, 
yet the best models do not fit the SED particularly well. 

Different dust compositions influence both the magni- 
tude and shape of the visibility spectrum. We now ex- 
plore dust made of cold and warm silicates by Ossenkopf 
et al. (1992) and amorphous carbon. Ossenkopf sili- 
cates have the strong advantage of a larger optical depth 
in the silicate feature than the DL grains. This prop- 
erty is advantageous in case of W33A, because it allows 
models with relatively moderate Ay to produce deep sili- 
cate absorption. Models with moderate extinction fit the 
shorter wavelength fluxes better. We find that the typical 
ISM ratio of 0.88 between graphite and silicates does not 
produce enough silicate absorption to fit W33A's strong 
absorption feature. A better correspondence is reached 
if this ratio is reduced to 0.50 (see also Churchwell et 
al. 1990). The warm Ossenkopf silicates bring the extra 
advantage of somewhat reducing the sublimation radius. 
Although the reduction in Ay, and a better correspon- 
dence to the silicate feature depth produces reasonable 
correspondence between the shallow density models with 
the observation, models that fit the data better are found 
by lowering the Teg of the star. The change of this pa- 
rameter produces a decrease of the size of the inner dust 
boundary, because of the lower dust opacity with tem- 
perature. We finally arrive at a simultaneous fit (Fig. [2]) 
to the red silicate wing, SED peak and visibilities for a 
central object with a relatively low Tcs of lO^'K (corre- 
sponding to a B9 supergiant, again consistent with the 
notion that an accreting star may be swollen) and a r~°'^ 
density law (Fig. [21). The dust sublimation radius for 
Td = 1500 K is found at 7 mas (26.5 AU). 

In summary, fitting ID DUSTY models to the visi- 



bilities and SED of W33A has shown that for nominal 
stellar parameters, size scales of the emission region are 
too large and produce the wrong trend with wavelength. 
Shallow radial density distributions produce this trend 
of size with wavelength and in order to fit the depth of 
the silicate feature as well, dust models with warm Os- 
senkopf silicates with an increased silicate over graphite 
ratio reproduce the observation best, provided the lumi- 
nosity and Teff of the central object are reduced. 

5. CONCLUSIONS 

We presented mid-IR high-resolution dispersed in- 
terferometric observations of the forming massive star 
W33A. The visibility spectrum indicates an equivalent 
Gaussian FWHM of the emitting region of ^ 120AU, at 
8^m, increasing to ~ 240AU, at 13/im. We interpreted 
the interferometric data with simple spherically symmet- 
ric DUSTY models representative of a (unresolved) star 
embedded in a dusty envelope, aiming for a simultane- 
ous fit to the SED, silicate profile and visibility spectrum. 
For any radial dust distribution we found that the canon- 
ical value of W33A's luminosity is not compatible with 
the visibilities. This is supported by MIPSGAL data. 
We found that even for a reduced luminosity, the model 
produces too large emitting regions causing the visibil- 
ities and SED to be mutually incompatible. Changing 
the dust composition improves the situation for an in- 
creased graphite to silicate ratio and using warm silicate 
optical constants from Ossenkopf et al. (1992). We re- 
sorted to a substantial lowering of the T^s in order to 
obtain a satisfactorily match between observables and 
models. Further coverage of the uv-plane will be very 
rewarding, constraining more appropriate, higher dimen- 
sionality models, which will eventually lead in a proper 
description of the circumstellar environment of an accret- 
ing massive star. A full 2D axi-symmetric treatment of 
this and other MIDI data will be the subject of a future 
paper. 
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